library(downloader) # pour télécharger des fichiers depuis un serveur
library(tidyverse) # pour simplifier la syntaxe
library(sf) # pour manipuler des objets SimpleFeature
library(sp) # pour manipuler des objets Spatial*DataFrame
library(spatstat) # pour calculer le lissage spatial
library(maptools) # pour convertir l'objet sp en objet ppp (spatstat)
library(cartography) # pour cartographier les résultats
library(cartogramR) # pour créer les cartogrammes
library(raster) # pour manipuler les rasters
library(terra) # pour manipuler les rasters
library(stars) #pour la vectorisation du raster
library(rmapshaper) # pour généraliser les contours (simplification des géométries)
library(dplyr) #pour manipuler des tableaux de données
library(openxlsx) #pour lire des fichiers EXCEL
library(ggplot2) #pour tracer des graphiques esthétiques
library(tidyr) #pour manipuler des matrices
library(rgdal) # pour manipuler des objets ayant une composante spatiale
library(rgeos) # pour calculer des centroides
library(spdep) # pour calculer l'auto-corrélation spatiale
library(geoR) # pour calculer le semi-variogramme empirique
library(magick) # pour créer un gif animé
library(av) # pour créer une vidéo mp4
library(knitr) # pour afficher les images exportées
library(filesstrings) # pour déplacer un fichier d'un dossier à un autre
download("https://github.com/ESO-Rennes/Animated-Cartograms/raw/main/Data_Pulsation.zip", destfile=paste0(getwd(),"/Data_Pulsation.zip"), mode="wb", overwrite = TRUE)
unzip(zipfile = "Data_Pulsation.zip", exdir = ".")
EMU2019 <- st_read(dsn = "Data_Pulsation", layer = "EMU2019", stringsAsFactors = FALSE)
## Reading layer `EMU2019' from data source
## `F:\Rennes\Recherche\ANR MoDurAL\Valorizacion\IJC\version anglaise\Derniers_traitements\Envoi Florent 12 01\Data_Pulsation'
## using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 56 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N
viajes <- read.xlsx("Data_Pulsation/ViajesEODH2019.xlsx") #134497 observations
t_col <- function(color, percent = 50, name = NULL) {
rgb.val <- col2rgb(color)
t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
max = 255,
alpha = (100 - percent) * 255 / 100,
names = name)
}
t_col_palette <- function(palette, percent = 50) {
palette_trans <- palette
for (j in 1:length(palette_trans)){
palette_trans[j] <- t_col(palette_trans[j], percent = trans)
}
return(palette_trans)
}
Calcul des durées de trajet (prise en compte des trajets qui commencent un jour et terminent le suivant)
viajes$duracion <- viajes$p31_hora_llegada - viajes$hora_inicio_viaje
for(i in 1:nrow(viajes)){
if(viajes$p31_hora_llegada[i]<viajes$hora_inicio_viaje[i]){
viajes$duracion[i] <- viajes$duracion[i]+1
}
}
Calcul du temps de trajet moyen en heures. Le temps de trajet moyen (aller-simple) est de 50 minutes
avg = 24*mean(viajes$duracion)
Tri pour écarter les trajets aberrants (0 minute)
viajes2 <- viajes %>% filter(duracion>0)
Traitement des trajets dont le point de départ ou d’arrivée se trouve en dehors de Bogotá DC. On crée une UTAM virtuelle regroupant les zones hors DC.
#Repérage des trajets qui sortent ou entrent dans le DC
excl <- c("N/A", "UPR1", "UPR2", "UPR3")
viajes2$utam_destino[is.na(viajes2$utam_destino)] <- "N/A"
viajes2$utam_origen[is.na(viajes2$utam_origen)] <- "N/A"
# création d'un UTAM virtuel "UTAM99"
for (i in 1:nrow(viajes2)){
if (viajes2$utam_destino[i] %in% excl){
viajes2$utam_destino[i] <- "UTAM999"
}
if (viajes2$utam_origen[i] %in% excl){
viajes2$utam_origen[i] <- "UTAM999"
}
}
On fait sortir la personne de l’UTAM d’origine à l’heure de départ et on la fait rentrer dans l’UTAM de destination à l’heure d’arrivée.
#Les heures et durées sont en fraction de la journée. On les multiplie par 24 pour les convertir en heures.
viajes2$h_arrivee <- 24*viajes2$p31_hora_llegada
viajes2$h_depart <- 24*viajes2$hora_inicio_viaje
#pas de temps de l'animation : défaut 15 minutes
pas_de_temps <- 0.25
#On arrondit les heures au quart d'heure précédent. Par exemple 8h27 --> 8h15
viajes2$h_arrivee_round <- pas_de_temps*(viajes2$h_arrivee-viajes2$h_arrivee%%pas_de_temps)%/%pas_de_temps
viajes2$h_depart_round <- pas_de_temps*(viajes2$h_depart-viajes2$h_depart%%pas_de_temps)%/%pas_de_temps
On filtre pour ne garder que les trajets qui font une boucle dans la journée afin d’éviter un décrochage lors de la lecture en boucle de l’animation. Pour chaque trajet d’un individu, on calcule la somme UTAM d’arrivée - UTAM de départ, puis on somme ces valeurs pour chaque individu. Les individus qui font une boucle sont repérés par la valeur 0. Les autres, que l’on souhaite exclure, ont une valeur strictement positive ou négative. Cette opération retire 5.7% des observations, ce que l’on considère acceptable.
viajes2$utam_origen_id <- as.numeric(gsub("UTAM", "", viajes2$utam_origen))
viajes2$utam_destino_id <- as.numeric(gsub("UTAM", "", viajes2$utam_destino))
viajes2$utam_test <- viajes2$utam_destino_id - viajes2$utam_origen_id
viajes_circulares_2 <- viajes2 %>%
group_by(id_hogar, id_persona) %>%
summarize(circular_2 = sum(utam_test))
viajes2 <- viajes2 %>%
left_join(viajes_circulares_2, by = c("id_hogar" = "id_hogar", "id_persona" = "id_persona"))
viajes2 <- viajes2 %>%
filter(circular_2 == 0)
On retire les trajets intra-UTAM qui sont considérés, à cette échelle, comme des trajets n’ayant pas eu lieu car ils ne contribuent pas à la déformation.
viajes2 <- viajes2 %>% filter(viajes2$utam_origen != viajes2$utam_destino)
On prépare un tableau contenant le solde de flux par UTAM
#solde de mouvements par utam pour une heure donnée pondérée par le facteur d'expansion
sorties_utam <- viajes2 %>%
group_by(h_depart_round,utam_origen) %>%
summarize(sorties = sum(f_exp))
entrees_utam <- viajes2 %>%
group_by(h_arrivee_round,utam_destino) %>%
summarize(entrees = sum(f_exp))
#création d'un tableau avec le solde de flux par UTAM
flux <- data.frame(EMU2019$UTAM)
names(flux)[names(flux) == 'EMU2019.UTAM'] <- 'UTAM'
flux <- flux %>%
full_join(sorties_utam, by = c("UTAM" = "utam_origen"))
flux <- flux %>%
full_join(entrees_utam, by = c("UTAM" = "utam_destino", "h_depart_round" = "h_arrivee_round"))
flux[is.na(flux)] <- 0
flux$solde <- flux$entrees - flux$sorties
names(flux)[names(flux) == 'h_depart_round'] <- 'h_round'
#tri chronologique par heure, puis par UTAM
flux <- arrange(flux, h_round, UTAM)
#conversion en matrice de flux
f <- pivot_wider(flux, id_cols = UTAM, names_from = h_round, values_from = solde)
f[is.na(f)] <- 0
Pour choisir l’heure de référence (où chacun est chez soi), on trace un histogramme de la répartition des heures de déplacement. On se rend compte que le creux de déplacement est maximal à 3h00. On va affecter la variable NUMPERSTOT, qui donne la population de nuit issue du recensement de 2018, à ce créneau.
#entrées
ggplot(flux, aes(x = h_round, y = entrees)) +
stat_summary(geom = "bar", position = "dodge", fun = sum)
#sorties
ggplot(flux, aes(x = h_round, y = sorties)) +
stat_summary(geom = "bar", position = "dodge", fun = sum)
#Création d'un tableau avec le stock de population par UTAM, heure par heure
stock <- f %>%
inner_join(EMU2019, by = 'UTAM', copy = TRUE) %>% #jointure avec EMU_2019
relocate('NUMPERSTOT', .after = '2.75') #positionnement de NUMPERSTOT au niveau de la tranche 2h45 du matin
#pour réordonner les colonnes en commençant à 3h00.
stock <- stock[,1:98]
stock1 <- stock[,1]
stock2 <- stock[,2:13]
stock3 <- stock[,14:98]
stock <-cbind(stock1, stock3, stock2)
#remplissage de la matrice de stock
for (i in 3:ncol(stock)){
stock[,i]<-stock[,i]+stock[,i-1]
}
#On remplace les éventuelles valeurs négatives par des zéros. #Aucune UTAM concernée a priori
for (i in 3:ncol(stock)){
for (j in 1:nrow(stock)){
if(stock[j,i]<0){
stock[j,i] <- 0
}
}
}
On transpose la table pour tracer des graphes de la population par UTAM au fil de la journée.
stock2 <- stock
rownames(stock2)<-stock2$UTAM #on renomme les lignes avec les codes des UTAM
stock2 <- stock2[,3:98] #on supprime les deux premières colonnes
tstock <- data.frame(t(stock2))
tstock$hours <- rownames(tstock)
tstock$rank <- c(1:96)
#à titre d'exemple, on trace le graphe pour deux UTAM : EL LUCERO et CHAPINERO
ggplot(tstock, aes(x = rank)) +
geom_line(aes(y=UTAM67, col="EL LUCERO")) +
geom_line(aes(y=UTAM99, col="CHAPINERO")) +
scale_colour_manual("", breaks = c("EL LUCERO", "CHAPINERO"), values = c("red", "blue")) +
scale_x_continuous(breaks=c(1+4*(0:23)), labels=c(3:23,0:2)) +
scale_y_continuous(breaks=c(50000,100000), labels=c("50 000","100 000")) +
labs(x="Time of day", y="Population") +
theme(legend.position = "left")
#Récupération des valeurs de stock dans le fond de carte
EMU2019_2 <- EMU2019[,c(2,55)] %>%
left_join(stock, by = "UTAM")
#on génère la liste des heures bien ordonnées
heures <- c(0.25*(12:95),0.25*(0:11))
#on génère des étiquettes à afficher sur les futures cartes (dans le titre) à partir de cette liste en les mettant au format anglophone
heures2 <- heures
for(i in 41:84){
heures2[i] <- heures2[i]-12
}
for(i in 85:88){
heures2[i] <- heures2[i]+12
}
heures_lab <- as.character(heures2)
for(i in 1:length(heures_lab)){
heures_lab[i] <- paste(as.character(heures2[i]%/%1),
".",as.character(60*heures2[i]%%1), sep = "")
if(i%%4 == 1){heures_lab[i] <- str_replace(heures_lab[i], fixed(".0"), fixed(".00"))}
}
for(i in 1:length(heures_lab)){
if(i %in% 37:84){
heures_lab[i] <- paste(heures_lab[i], "pm", sep = " ")
}
else{
heures_lab[i] <- paste(heures_lab[i], "am", sep = " ")
}
}
#on renomme les colonnes des heures
for (i in 1:ncol(EMU2019_2)){
names(EMU2019_2)[names(EMU2019_2) == heures[i]] <- paste("stock",heures[i], sep = "_")
}
Tests d’export de cartogrammes sur des heures individuelles, dans le cas présent à 3h00. On fixe une erreur absolue d’une valeur de 10 hectares, soit 100 000 m². On limite le calcul à 100 itérations.
EMU2019_carto <- cartogramR(EMU2019_2, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 100, absrel = FALSE, abserror = 100000, L = 256))
titre <- paste("Population per UTAM in Bogotá", heures_lab[1], sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(EMU2019_carto, color = "black")
title(main = titre, cex.main = 2)
dev.off()
## png
## 2
include_graphics(paste(titre,".png"))
#conversion en fichier sf
EMU2019_carto_sf <- as.sf(EMU2019_carto)
On va calculer les facteurs multiplicatifs à partir de la table de stock. On remplit le tableau colonne par colonne en divisant le stock pour chaque tranche horaire par le stock de référence NUMPERSTOT.
variation <- stock
for (i in 3:ncol(variation)){
variation[,i]<-variation[,i]/variation[,2]
#on renomme les colonnes pour les distinguer des variables de stock lors de la jointure
names(variation)[names(variation) == heures[i]] <- paste("var",heures[i], sep = "_")
}
names(variation)[names(variation) == heures[1]] <- paste("var",heures[1], sep = "_")
names(variation)[names(variation) == heures[2]] <- paste("var",heures[2], sep = "_")
#Jointure avec le fichier des géométries
EMU2019_3 <- EMU2019_2 %>%
left_join(variation, by = "UTAM")
Choix des seuils pour 6 classes
bks <- c(0,0.7,0.9,1.1,2,5,50)
Création d’une palette de couleurs (double gradation harmonique)
cols <- c("#303b98","#00beed","#fafaaa","#f2bf26","#e96831","#bb2b30")
#facteur de transparence
trans <- 50
#application du facteur de transparence
for (j in 1:6){
cols[j] <- t_col(cols[j], percent = trans)
}
titre <- paste("Multiplication factor with respect to night population", heures_lab[37], sep = " - ")
png(paste(titre, ".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3))
choroLayer(x = EMU2019_3, var = paste("var", heures[37], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1.5, legend.title.cex = 2)
title(main = paste("Multiplication factor with respect to night population", heures_lab[37], sep = "\n"), cex.main = 2)
dev.off()
## png
## 2
include_graphics(paste(titre,".png"))
#Pour convertir la couche des UTAM au format spatial object (format requis par le package geoR)
UTAM <- as(EMU2019_3, "Spatial")
#Pour calculer les centroides des UTAM (indispensable pour le calcul sur les plus proches voisins)
UTAMCentroids <- gCentroid(UTAM,byid=TRUE)
#Pour récupérer les données initiales des UTAM sur les centroides
UTAMCentroids <- SpatialPointsDataFrame(UTAMCentroids, UTAM@data)
#Calcul sur les plus proches voisins
listPPV <- knearneigh(UTAMCentroids@coords, k = 1) # pour connaître le plus proche voisin de chaque UTAM
PPV <- knn2nb(listPPV, row.names = UTAM$UTAM) # pour convertir l'objet knn en objet nb
distPPV <- nbdists(PPV, UTAMCentroids@coords) # pour connaître la distance entre plus proches voisins
print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 553.5302 1365.71 1587.616 2181.223 1854.208 10899.47
hist(unlist(distPPV), breaks = 20,
main = "Distance to the closest neighbour",
col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")
La plupart des UTAM ont au moins un voisin dans un rayon de 1500m.
#pour convertir les UTAM en objects nb
nbUTAM <- poly2nb(pl = UTAM,
row.names = UTAM$UTAM,
snap = 50,
queen = TRUE)
# pour identifier les UTAM sans voisins topologiques
summary(nbUTAM)
## Neighbour list object:
## Number of regions: 132
## Number of nonzero links: 632
## Percentage nonzero weights: 3.627181
## Average number of links: 4.787879
## 16 regions with no links:
## UTAM89 UTAM563 UTAM540 UTAM580 UTAM640 UTAM650 UTAM660 UTAM670 UTAM680
## UTAM690 UTAM700 UTAM600 UTAM630 UTAM620 UTAM610 UTAM590
## Link number distribution:
##
## 0 1 2 3 4 5 6 7 8 9 10
## 16 3 4 12 19 19 25 16 11 6 1
## 3 least connected regions:
## UTAM68 UTAM500 UTAM520 with 1 link
## 1 most connected region:
## UTAM112 with 10 links
# création d'une liste des 16 UTAM sans voisins fournis par l'instruction précédente
UTAM_isolees <- c("UTAM89",
"UTAM563",
"UTAM540",
"UTAM580",
"UTAM640",
"UTAM650",
"UTAM660",
"UTAM670",
"UTAM680",
"UTAM690",
"UTAM700",
"UTAM600",
"UTAM630",
"UTAM620",
"UTAM610",
"UTAM590")
#suppression des "UTAM_isolees" sans quoi le calcul de l'indice de Moran n'est pas possible
UTAM <- UTAM[which(!UTAM$UTAM %in% UTAM_isolees),]
EMU2019_3 <- EMU2019_3[which(!EMU2019_3$UTAM %in% UTAM_isolees),]
nbUTAM <- poly2nb(pl = UTAM,
row.names = UTAM$UTAM,
snap = 50,
queen = TRUE)
#Calcul du test de Moran pour chaque heure
moran <- c(0*(1:96))
for(i in 1:96){
m <- moran.test(UTAM@data[,100+i], listw = nb2listw(nbUTAM)) #dans la table UTAM, les colonnes contenant les facteurs multiplicatifs sont les 101 à 196
moran[i] <- as.numeric(m$estimate[1])
}
#affichage
df <- data.frame(cbind(moran,heures))
ggplot(df) +
theme_classic() +
geom_point(aes(x = heures, y=moran), shape = 21, fill = "white", color = "black", size = 2, stroke = 0.5) +
labs(x="Time of day", y="Moran") +
theme(legend.position = "left",
panel.border = element_rect(colour = "black", fill=NA, size=0.5))
### Variogramme
#pour calculer les pseudo-centroides des UTAM sans les UTAM isolées (indispensable pour le calcul du semi-variogramme)
UTAMCentroids <- gPointOnSurface(UTAM,byid=TRUE)
#pour récupérer les données initiales des communes sur les centroides
UTAMCentroids <- SpatialPointsDataFrame(UTAMCentroids, UTAM@data)
#pour convertir le SpatialPointsDataFrame en objet geodata
UTAMCentroids.geodata <- as.geodata(UTAMCentroids, data.col = "var_12")
#pour calculer le semi-variogramme empirique
vario.ex<- variog(UTAMCentroids.geodata, bin.cloud=TRUE, option = "bin")
## variog: computing omnidirectional variogram
plot(vario.ex, main = "Semivariogram of the multiplication factor at 12.00 pm in function of the distance", cex.main = 1)
lines(vario.ex, type ="l", lty = 2, col="red")
### Carte lissée à 12h00
#pour définir le contour de Bogotá comme emprise pour le lissage
Emprise <- as.owin(gBuffer(UTAM, width=0))
#pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM$var_12)
#pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha) --> calcul long
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
#Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(0, 0, 0, 0))
#Affichage de toutes les UTAM en arrière plan en centrant la carte sur Bogotá.
plot(st_geometry(EMU2019_3), col = NA, border = "black", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
#Affichage de la surface lissée
plot(cartelissee.raster, breaks = bks, col=cols, add = T, legend=F)
#Affichage des limites des UTAM
plot(st_geometry(EMU2019_3), border = "black", lwd = 0.05, lty=3, add = T)
legendChoro(
pos = "bottomleft",
title.txt = "Multiplication factor",
breaks = bks,
nodata = FALSE,
values.rnd = 2,
col = cols,
cex = 1.2,
values.cex = 1,
title.cex = 1
)
title("Population multiplication factor at 12.00 pm \n Smoothing radius 1000m", cex.main = 1, line = -2)
Le choix du rayon de lissage est important. Après plusieurs essais, on retient un rayon de 1000m qui est un bon compromis entre la prise en compte de suffisamment de voisins pour une UTAM donnée et la représentation la plus fidèle de la dynamique des différentes UTAM. Ci-dessous, la comparaison entre un rayon de lissage de 1000m et un rayon de 1500m.
##Isolement de l’UTAM “El Rincon de Suba” pour l’afficher en tant que légende sur le côté de la carte et représenter l’échelle de la déformation
EMU2019_carto <- cartogramR(EMU2019_2, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 100, absrel = FALSE, abserror = 100000, L = 256))
EMU2019_carto_sf <- as.sf(EMU2019_carto)
RINCON <- EMU2019_carto_sf[which(EMU2019_carto_sf$UTAM == "UTAM28"),]
st_geometry(RINCON)[[1]] = st_geometry(RINCON)[[1]] + st_point(c(-25000,-10000))
RINCON$legend <- paste(RINCON$NUMPERSTOT, "people", sep = "\n")
Pour ce faire, on réalise une interpolation entre la carte de Bogotá non déformée et la carte déformée selon la population de nuit à 3h00.
#Déformation du fond de carte initial pour le cartogramme seul
dir.create("Animated_Cartogram")
#Calcul des surfaces
EMU2019_area <- EMU2019[which(!EMU2019$UTAM %in% UTAM_isolees),]
EMU2019_area$Area <- as.numeric(st_area(st_geometry(EMU2019_area)))
#calcul de la population totale
Population <- sum(EMU2019_area$NUMPERSTOT)
#calcul de la superficie totale
Superficie <- sum(EMU2019_area$Area)
#nombre d'images intermédiaires pour la déformation
imgs <- 8
#Ventilation de la population totale proportionnellement à la superficie des UTAM pour l'image 0
EMU2019_area$NUMPERSTOT_0 <- EMU2019_area$Area*Population/Superficie
#Champ que l'on va faire varier
EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT
#création des images intermédiaires par une interpolation linéaire (plus lisse qu'en utilisant le TCAM, qui donnerait une croissance trop rapide sur la fin et un arrêt brutal)
Emprise <- as.owin(gBuffer(UTAM, width=0))
par(mar = c(1, 1, 1, 1))
#affichage de la carte non déformée
png(paste0(getwd(),"/Animated_Cartogram/", "deformation 0.png"), width = 960, height = 960)
plot(st_geometry(EMU2019_area), border = "black", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
dev.off()
## png
## 2
for (i in 1:imgs){
EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT*i/imgs + EMU2019_area$NUMPERSTOT_0*(imgs-i)/imgs
carto <- cartogramR(EMU2019_area, count = "NUMPERSTOT_i", options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
png(paste0(getwd(),"/Animated_Cartogram/", paste0("deformation ", i,".png")), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = NA, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
carto_sf <- as.sf(carto)
plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
dev.off()
}
#Déformation du fond de carte initial pour le cartogramme seconde variable.
dir.create("Choropleth_animated_cartogram")
#Calcul des surfaces
EMU2019_area <- EMU2019[which(!EMU2019$UTAM %in% UTAM_isolees),]
EMU2019_area$Area <- as.numeric(st_area(st_geometry(EMU2019_area)))
#calcul de la population totale
Population <- sum(EMU2019_area$NUMPERSTOT)
#calcul de la superficie totale
Superficie <- sum(EMU2019_area$Area)
#nombre d'images intermédiaires pour la déformation
imgs <- 8
#Ventilation de la population totale proportionnellement à la superficie des UTAM pour l'image 0
EMU2019_area$NUMPERSTOT_0 <- EMU2019_area$Area*Population/Superficie
#Champ que l'on va faire varier
EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT
#création des images intermédiaires par une interpolation linéaire (plus lisse qu'en utilisant le TCAM, qui donnerait une croissance trop rapide sur la fin et un arrêt brutal)
Emprise <- as.owin(gBuffer(UTAM, width=0))
par(mar = c(1, 1, 1, 1))
#affichage de la carte non déformée
png(paste0(getwd(),"/Choropleth_animated_cartogram/", "deformation 0.png"), width = 960, height = 960)
plot(st_geometry(EMU2019_area), col = NA, border = "black", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
dev.off()
## png
## 2
for (i in 1:imgs){
EMU2019_area$NUMPERSTOT_i <- EMU2019_area$NUMPERSTOT*i/imgs + EMU2019_area$NUMPERSTOT_0*(imgs-i)/imgs
carto <- cartogramR(EMU2019_area, count = "NUMPERSTOT_i", options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
png(paste0(getwd(),"/Choropleth_animated_cartogram/", paste0("deformation ", i,".png")), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
carto_sf <- as.sf(carto)
plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
dev.off()
}
# duplicate first initial images
dir.create("Smoothed_binned_colored_animated_cartogram")
initial_files <- list.files(paste0(getwd(),"/Choropleth_animated_cartogram/"))
file.copy(from = paste0(getwd(),"/Choropleth_animated_cartogram/", initial_files),
to = paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram/", initial_files))
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
dir.create("Smoothed_continuous_colored_animated_cartogram")
file.copy(from = paste0(getwd(),"/Choropleth_animated_cartogram/", initial_files),
to = paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram/", initial_files))
## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
setwd(paste0(getwd(),"/Animated_Cartogram"))
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
#niveau de transparence (%)
for (i in c(90, 80, 70, 60, 50)) {
trans <- i
Emprise <- as.owin(gBuffer(UTAM, width=0))
par(mar = c(1, 1, 1, 1))
# Affichage de toutes les UTAM en arrière-plan
titre <- paste("transition", 100 - trans, "percent", sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "white", bg = "white", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
par(col.main = paste("grey", 1.8*i-90))
title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2)
# conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
#Affichage de l'UTAM EL RINCON dans un coin
plot(st_geometry(RINCON), col = NA, border = paste("grey", 1.8*i-90), add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0, col = paste("grey", 1.8*i-90))
#Affichage de la carte non déformée
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = paste("grey", 0.6*i+30), add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, col = paste("grey", 1.8*i-90), outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))
Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
EMU2019_carto_sf <- as.sf(carto)
for (i in c(90,80,70,60,50)) {
trans <- i
#Affichage du fond
par(mar = c(1,1,1,1))
titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
choroLayer(x = EMU2019_carto_sf, var = paste("var", heures[1], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1, legend.title.cex = 1.5, add = TRUE)
title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2)
#Affichage de l'UTAM EL RINCON dans un coin
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#Affichage de la carte non déformée
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(n = 1)
background <- t_col("grey90",200-2*i)
mat <- matrix(c(0,1.0,0,1), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), col = NA, border = background, bg = background, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
plot(st_geometry(EMU2019_carto_sf), border = "black", lwd = 1, add = TRUE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,196])
# pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha) --> calcul long
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
# Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
# reclassification de la surface lissée préalablement calculée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
#vectorisation de la surface reclassée
cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
# pour trier les enregistrements de la colonne layer
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
# création d'une liste de classes présentes à l'heure donnée
liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
# Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
# Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
# pour généraliser les contours de cartelissee.vecteur.drapee (on garde 5% des sommets initiaux)
cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.05, keep_shapes = TRUE)
#niveau de transparence (%)
for (i in c(90, 80, 70, 60, 50)) {
trans <- i
cols_2 <- cols
for (j in 1:10){
cols_2[j] <- t_col(cols[j], percent = trans)
}
Emprise <- as.owin(gBuffer(UTAM, width=0))
par(mar = c(1, 1, 1, 1))
# création d'un tableau de correspondances avec pour chaque classe son code couleur
Correspondance_palette <- as.data.frame(cbind(c(1:6), cols_2))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
# création d'une palette de couleurs correspondant uniquement aux classes présentes
couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
# Affichage de toutes les UTAM en arrière-plan
titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
# Affichage de la surface lissée déformée
typoLayer(
x = cartelissee.vecteur.drapee,
var="layer",
col = as.character(couleurs_classes_presentes$cols),
lwd = 0.1,
border = as.character(couleurs_classes_presentes$cols),
legend.pos = "n",
add=T)
legendChoro(
pos = "bottomleft",
title.txt = "Multiplication factor",
breaks = c("","0.72","0.9","1.1","2","5",""),
nodata = FALSE,
values.rnd = 2,
col = cols,
cex = 1.2,
values.cex = 1,
title.cex = 1.5
)
par(col.main = paste("grey", 1.8*i-90))
title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2)
# conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
# Affichage
plot(st_geometry(carto_sf), border = paste("grey", 200-2*i), lwd = 1, add = TRUE)
#Affichage des contours des localités
Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
plot(st_geometry(Localidades_carto_2), col = NA, border = "black", lwd = 3.25 - i/40, add = TRUE)
#Affichage de l'UTAM EL RINCON en légende
plot(st_geometry(RINCON), col = NA, border = paste("grey", 1.8*i-90), add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0, col = paste("grey", 1.8*i-90))
#Affichage de la carte non déformée dans un carton en bas à droite
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = paste("grey", 0.6*i+30), add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, col = paste("grey", 1.8*i-90), outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
# paramétrage de la discrétisation (définition des seuils ou bornes de classes)
bks_pseudo_continuous <- c(seq(from=0, to=3, by = 0.025), seq(from=3.5, to=50, by = 0.5))
n0<-length(seq(from=0, to=0.6, by = 0.025))
n1<-length(seq(from=0.625, to=0.8, by = 0.025))
n2<-length(seq(from=0.825, to=0.975, by = 0.025))
n3<-length(seq(from=1, to=1.5, by = 0.025))
n4<-length(seq(from=1.525, to=3, by = 0.025))
n5<-length(seq(from=3.5, to=10, by = 0.5))
n6<-length(seq(from=10.5, to=50, by = 0.5))
#palette de couleur (double gradation harmonique) et facteur de transparence
cols_continuous <- c("#1d235c", "#303b98","#00beed","#fafaaa","#f2bf26","#e96831","#bb2b30","#7f0000")
trans <- 50
pal0 <- colorRampPalette(colors = cols_continuous[1:2], interpolate = "linear")
palette0 <- pal0(n0)
palette0 <- t_col_palette(palette0, trans)
pal1 <- colorRampPalette(colors = cols_continuous[2:3], interpolate = "linear")
palette1 <- pal1(n1)
palette1 <- t_col_palette(palette1, trans)
pal2 <- colorRampPalette(colors = cols_continuous[3:4], interpolate = "linear")
palette2 <- pal2(n2)
palette2 <- t_col_palette(palette2, trans)
pal3 <- colorRampPalette(colors = cols_continuous[4:5], interpolate = "linear")
palette3 <- pal3(n3)
palette3 <- t_col_palette(palette3, trans)
pal4 <- colorRampPalette(colors = cols_continuous[5:6], interpolate = "linear")
palette4 <- pal4(n4)
palette4 <- t_col_palette(palette4, trans)
pal5 <- colorRampPalette(colors = cols_continuous[6:7], interpolate = "linear")
palette5 <- pal5(n5)
palette5 <- t_col_palette(palette5, trans)
pal6 <- colorRampPalette(colors = cols_continuous[7:8], interpolate = "linear")
palette6 <- pal6(n6)
palette6 <- t_col_palette(palette6, trans)
palette <- c(palette0, palette1, palette2, palette3, palette4, palette5, palette6)
#preparation de la légende continue
leg1 <- pal1(10)
leg1 <- t_col_palette(leg1, trans)
leg2 <- pal2(10)
leg2 <- t_col_palette(leg2, trans)
leg3 <- pal3(10)
leg3 <- t_col_palette(leg3, trans)
leg4 <- pal4(10)
leg4 <- t_col_palette(leg4, trans)
leg5 <- pal5(10)
leg5 <- t_col_palette(leg5, trans)
leg <- cbind(leg1, leg2, leg3, leg4, leg5)
Emprise <- as.owin(gBuffer(UTAM, width=0))
Correspondance_palette <- as.data.frame(cbind(c(1:(n0+n1+n2+n3+n4+n5+n6)), palette))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
names(Correspondance_palette)[names(Correspondance_palette) == "V2"] <- "palette"
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[1], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
#pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,196])
#pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
# pour découper le raster sur l'emprise des UTAM
cartelissee.raster <- mask(cartelissee.raster, UTAM)
#reclassification de la surface lissée préalablement calculée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks_pseudo_continuous)
#vectorisation de la surface reclassée
cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
#pour trier les enregistrements de la colonne layer
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
#création d'une liste de classes présentes à l'heure donnée
liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
l0 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer <= n0]
l1 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+1):(n0+n1)]
l2 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+1):(n0+n1+n2)]
l3 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+1):(n0+n1+n2+n3)]
l4 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+1):(n0+n1+n2+n3+n4)]
l5 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+n4+1):(n0+n1+n2+n3+n4+n5)]
l6 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer > (n0+n1+n2+n3+n4+n5)]
#création d'une palette de couleurs correspondant uniquement aux classes présentes
couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
# Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
#Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
#pour généraliser les contours de cartelissee.vecteur.drapee (on garde 10% des sommets initiaux)
cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.1, keep_shapes = TRUE)
#niveau de transparence (%)
for (i in c(50,60,70,80,90)) {
trans <- i
#Affichage du fond
par(mar = c(1,1,1,1))
titre <- paste("transition", 100 - trans, "pourcent", sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
#Affichage de la surface lissée déformée (valeurs jusqu'à 0.6)
if(length(l0)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l0, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l0]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeurs de 0.6 à 0.8)
if(length(l1)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l1, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l1]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 0.8 à 1)
if(length(l2)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l2, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l2]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 1 à 1.5)
if(length(l3)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l3, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l3]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 1.5 à 3)
if(length(l4)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l4, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l4]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 3 à 10)
if(length(l5)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l5, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l5]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur à partir de 10)
if(length(l6)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l6, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l6]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#seulement pour l'échelle de la légende, sans les couleurs
legendChoro(
pos = "bottomleft",
title.txt = "Multiplication factor",
breaks = c("","0.72","0.9","1.1","2","5",""),
nodata = FALSE,
values.rnd = 2,
col = c("grey90","grey90","grey90","grey90","grey90","grey90"),
cex = 1.2,
values.cex = 1,
title.cex = 1.5,
border = "grey90"
)
title(main = "Population per UTAM in Bogotá - 3.00 am", cex.main = 2, line = -2)
# conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
#Affichage des contours des localités
Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
#Affichage de l'UTAM EL RINCON dans un coin
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#affichage de la légende continue
mat <- matrix(c(0.013,0.14,0.013,0.297), nrow=1)
split.screen(figs = mat, erase=FALSE)
image(1, 1:length(leg), t(as.matrix(1:length(leg))), col = leg, xlab = "", ylab = "", xaxt = "n",
yaxt = "n",bty = "n")
close.screen(n = 1)
#Affichage de la carte non déformée
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(n = 1)
background <- t_col("grey90",200-2*i)
mat <- matrix(c(0,1.0,0,1), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), col = NA, border = background, bg = background, xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
plot(st_geometry(carto_sf), border = paste("grey", 200-2*i), lwd = 1, add = TRUE)
plot(st_geometry(Localidades_carto_2), col = NA, border = "black", lwd = 3.25 - i/40, add = TRUE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Animated_Cartogram"))
Emprise <- as.owin(gBuffer(UTAM, width=0))
for (i in 1:96){
par(mar = c(1, 1, 1, 1))
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
# Affichage du fond
plot(st_geometry(EMU2019_3), col = NA, border = "white", bg = "white", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
title(main = titre, cex.main = 2, line = -2)
# conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
plot(st_geometry(carto_sf), border = "black", lwd = 1, add = TRUE)
#Affichage de l'UTAM EL RINCON en légende
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#Affichage de la carte non déformée dans un carton en bas à droite
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))
Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
for (i in 1:96){
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
EMU2019_carto_sf <- as.sf(carto)
par(mar = c(1, 1, 1, 1))
titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
choroLayer(x = EMU2019_carto_sf, var = paste("var", heures[i], sep = "_"), breaks = bks, col = cols, legend.title.txt = "Multiplication factor", legend.values.rnd = 2, legend.values.cex = 1, legend.title.cex = 1.5, add = TRUE)
title(main = titre, cex.main = 2, line = -2)
#Affichage de l'UTAM EL RINCON dans un coin
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#Affichage de la carte non déformée
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
Emprise <- as.owin(gBuffer(UTAM, width=0))
Correspondance_palette <- as.data.frame(cbind(c(1:6), cols))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
for (i in 1:96){
par(mar = c(1, 1, 1, 1))
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
#pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,100+i])
#pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
#reclassification de la surface lissée préalablement calculée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks)
#vectorisation de la surface reclassée
cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE)
#pour trier les enregistrements de la colonne layer
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
#création d'une liste de classes présentes à l'heure donnée
liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
#création d'une palette de couleurs correspondant uniquement aux classes présentes
couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
# Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
#Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
#pour généraliser les contours de cartelissee.vecteur.drapee (on garde 5% des sommets initiaux)
cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.05, keep_shapes = TRUE)
#Affichage du fond
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
#Affichage de la surface lissée déformée
typoLayer(
x = cartelissee.vecteur.drapee,
var="layer",
col = as.character(couleurs_classes_presentes$cols),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
legendChoro(
pos = "bottomleft",
title.txt = "Multiplication factor",
breaks = c("","0.72","0.9","1.1","2","5",""),
nodata = FALSE,
values.rnd = 2,
col = cols,
cex = 1.2,
values.cex = 1,
title.cex = 1.5
)
title(main = titre, cex.main = 2, line = -2)
#conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
#Affichage des limites des UTAM déformées
plot(st_geometry(carto_sf), border = "white", lwd = 1, add = TRUE)
#Affichage des zones d'étude MODURAL (optionnel)
#plot(st_geometry(carto_sf_modural), border = "black", lwd = 3, add = TRUE)
#plot(st_geometry(carto_sf_modural), border = "white", lwd = 1, add = TRUE)
#Affichage des contours des localités
Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
plot(st_geometry(Localidades_carto_2), col = NA, lwd = 2, add = TRUE)
#Affichage de l'UTAM EL RINCON en légende
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#Affichage de la carte non déformée dans un carton en bas à droite
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
Emprise <- as.owin(gBuffer(UTAM, width=0))
Correspondance_palette <- as.data.frame(cbind(c(1:(n0+n1+n2+n3+n4+n5+n6)), palette))
names(Correspondance_palette)[names(Correspondance_palette) == "V1"] <- "Id_Classe"
names(Correspondance_palette)[names(Correspondance_palette) == "V2"] <- "palette"
for (i in 1:96){
par(mar = c(1, 1, 1, 1))
carto <- cartogramR(EMU2019_3, count = paste("stock",heures[i], sep = "_"), options = list(maxit = 20, absrel = FALSE, abserror = 100000, L = 256))
titre <- paste("Population per UTAM in Bogotá", heures_lab[i], sep = " - ")
png(paste(titre,".png"), width = 960, height = 960)
#pour creer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (facteur multiplicatif à 12H)
UTAM.ppp <- ppp(UTAMCentroids@coords[,1], UTAMCentroids@coords[,2], window = Emprise, marks = UTAM@data[,100+i])
#pour calculer la surface lissée (rayon lissage : 1 km et résolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(UTAM.ppp, kernel = "gaussian", sigma = 1000, weights = UTAM.ppp$marks, eps=c(100,100))
#Conversion de la surface lissée au format raster
cartelissee.raster <- raster(cartelissee)
crs(cartelissee.raster) <- st_crs(UTAM)$srid # pour spécifier un SCR à l'objet raster
# pour découper le raster sur l'emprise des UTAM
cartelissee.raster <- mask(cartelissee.raster, UTAM)
#reclassification de la surface lissée préalablement calculée
cartelissee.reclass <- cut(cartelissee.raster, breaks = bks_pseudo_continuous)
#vectorisation de la surface reclassée
cartelissee.vecteur <- sf::st_as_sf(stars::st_as_stars(cartelissee.reclass), as_points = FALSE, merge = TRUE) # requires the sf, sp, raster and stars packages
#pour trier les enregistrements de la colonne layer
cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$layer, decreasing = FALSE), ]
#création d'une liste de classes présentes à l'heure donnée
liste_classes_presentes <- as.character(cartelissee.vecteur$layer)
l0 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer <= n0]
l1 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+1):(n0+n1)]
l2 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+1):(n0+n1+n2)]
l3 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+1):(n0+n1+n2+n3)]
l4 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+1):(n0+n1+n2+n3+n4)]
l5 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer %in% (n0+n1+n2+n3+n4+1):(n0+n1+n2+n3+n4+n5)]
l6 <- cartelissee.vecteur$layer[cartelissee.vecteur$layer > (n0+n1+n2+n3+n4+n5)]
#création d'une palette de couleurs correspondant uniquement aux classes présentes
couleurs_classes_presentes <- Correspondance_palette[which(Correspondance_palette$Id_Classe %in% liste_classes_presentes),]
# Pour rogner légèrement cartelissee.vecteur (-50 mètres) pour que l'emprise de cet objet soit comprise à l'intérieur de celle du cartogramme (opération à n'effectuer que si la commande suivante ne fonctionne pas)
cartelissee.vecteur <- st_crop(cartelissee.vecteur, ext(UTAM)-50)
#Pour déformer la surface lissée suivant les mêmes paramètres que ceux ayant servi à construire le cartogramme
cartelissee.vecteur.drapee <- geom_cartogramR(cartelissee.vecteur, carto)
#pour généraliser les contours de cartelissee.vecteur.drapee (on garde 10% des sommets initiaux)
cartelissee.vecteur.drapee <- ms_simplify(cartelissee.vecteur.drapee, keep = 0.1, keep_shapes = TRUE)
#Affichage du fond
plot(st_geometry(EMU2019_3), col = NA, border = "grey90", bg = "grey90", xlim = c(st_bbox(UTAM)[1], st_bbox(UTAM)[3]), ylim = c(st_bbox(UTAM)[2], st_bbox(UTAM)[4]))
#Affichage de la surface lissée déformée (valeurs jusqu'à 0.6)
if(length(l0)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l0, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l0]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeurs de 0.6 à 0.8)
if(length(l1)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l1, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l1]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 0.8 à 1)
if(length(l2)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l2, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l2]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 1 à 1.5)
if(length(l3)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l3, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l3]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 1.5 à 3)
if(length(l4)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l4, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l4]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur de 3 à 10)
if(length(l5)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l5, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l5]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#Affichage de la surface lissée déformée (valeur à partir de 10)
if(length(l6)>0){
typoLayer(
x = cartelissee.vecteur.drapee[cartelissee.vecteur.drapee$layer %in% l6, ],
var="layer",
col = as.character(couleurs_classes_presentes$palette[couleurs_classes_presentes$Id_Classe %in% l6]),
lwd = 0.1,
border = NA,
legend.pos = "n",
add=T)
}
#seulement pour l'échelle de la légende, sans les couleurs
legendChoro(
pos = "bottomleft",
title.txt = "Multiplication factor",
breaks = c("","0.72","0.9","1.1","2","5",""),
nodata = FALSE,
values.rnd = 2,
col = c("grey90","grey90","grey90","grey90","grey90","grey90"),
cex = 1.2,
values.cex = 1,
title.cex = 1.5,
border = "grey90"
)
title(main = titre, cex.main = 2, line = -2)
# conversion de l'objet cartogram en objet sf
carto_sf <- as.sf(carto)
# Affichage
plot(st_geometry(carto_sf), border = "white", lwd = 1, add = TRUE)
#Affichage des contours des localités
Localidades_2 <- EMU2019_3 %>% group_by(LocMuni) %>% summarize()
Localidades_carto_2 <- geom_cartogramR(Localidades_2, carto)
plot(st_geometry(Localidades_carto_2), col = NA, lwd = 2, add = TRUE)
#Affichage de l'UTAM EL RINCON dans un coin
plot(st_geometry(RINCON), col = NA, border = "black", add = TRUE)
labelLayer(RINCON, txt = "legend", halo = TRUE, bg = "white", r = 0.05, cex = 1.2, pos = 3, font = 2, offset = 0)
#affichage de la légende continue
mat <- matrix(c(0.013,0.14,0.013,0.297), nrow=1)
split.screen(figs = mat, erase=FALSE)
image(1, 1:length(leg), t(as.matrix(1:length(leg))), col = leg, xlab = "", ylab = "", xaxt = "n",
yaxt = "n",bty = "n")
close.screen(n = 1)
#Affichage de la carte non déformée
mat <- matrix(c(0.7,1.0,0.1,0.4), nrow=1)
split.screen(figs = mat, erase=FALSE)
plot(st_geometry(EMU2019_3), border = "grey60", add = TRUE)
mtext("Initial undistorted map", side = 1, line = 3, outer = FALSE)
close.screen(all = TRUE)
dev.off()
}
dir.create("Animated_Cartogram_Export")
#chargement des informations des fichiers à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Animated_Cartogram"), full.names = TRUE))
#tri par date de création
details = details[with(details, order(as.POSIXct(ctime))), ]
list <- rep("0",110)
for (i in 1:10){
list[i] <- paste("00",i, sep = "")
}
for (i in 10:110){
if (i>=10 & i<100){
list[i] <- paste("0",i, sep = "")
} else {
list[i] <- paste(i)
}
}
file.rename(rownames(details), paste0(getwd(),"/Animated_Cartogram/anim", list[1:110],".png"))
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE
# rechargement des informations des fichiers renommés à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Animated_Cartogram"), full.names = TRUE))
details <- details[15:110,]
imgs = rownames(details)
img_list <- lapply(imgs, image_read)
img_joined <- image_join(img_list)
img_animated <- image_animate(img_joined, fps = 20)
#export en gif
image_write(image = img_animated, path = "Animated_Cartogram_Export/animation.gif")
include_graphics(paste0(getwd(),"/Animated_Cartogram_Export/animation.gif"))
#export en mp4
rootwd <- getwd()
setwd(paste0(getwd(),"/Animated_Cartogram"))
av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")
## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Animated_Cartogram\\animation.mp4"
file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Animated_Cartogram_Export"))
dir.create("Choropleth_animated_cartogram_Export")
#chargement des informations des fichiers à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Choropleth_animated_cartogram"), full.names = TRUE))
#tri par date de création
details = details[with(details, order(as.POSIXct(ctime))), ]
list <- rep("0",110)
for (i in 1:10){
list[i] <- paste("00",i, sep = "")
}
for (i in 10:110){
if (i>=10 & i<100){
list[i] <- paste("0",i, sep = "")
} else {
list[i] <- paste(i)
}
}
file.rename(rownames(details), paste0(getwd(),"/Choropleth_animated_cartogram/anim", list[1:110],".png"))
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE
# rechargement des informations des fichiers renommés à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Choropleth_animated_cartogram"), full.names = TRUE))
details <- details[15:110,]
imgs = rownames(details)
img_list <- lapply(imgs, image_read)
img_joined <- image_join(img_list)
img_animated <- image_animate(img_joined, fps = 20)
#export en gif
image_write(image = img_animated, path = "Choropleth_animated_cartogram_Export/animation.gif")
include_graphics(paste0(getwd(),"/Choropleth_animated_cartogram_Export/animation.gif"))
#export en mp4
rootwd <- getwd()
setwd(paste0(getwd(),"/Choropleth_animated_cartogram"))
av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")
## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Choropleth_animated_cartogram\\animation.mp4"
file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Choropleth_animated_cartogram_Export"))
dir.create("Smoothed_binned_colored_animated_cartogram_Export")
#chargement des informations des fichiers à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"), full.names = TRUE))
#tri par date de création
details = details[with(details, order(as.POSIXct(ctime))), ]
list <- rep("0",110)
for (i in 1:10){
list[i] <- paste("00",i, sep = "")
}
for (i in 10:110){
if (i>=10 & i<100){
list[i] <- paste("0",i, sep = "")
} else {
list[i] <- paste(i)
}
}
file.rename(rownames(details), paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram/anim", list[1:110],".png"))
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE
# rechargement des informations des fichiers renommés à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"), full.names = TRUE))
details <- details[15:110,]
imgs = rownames(details)
img_list <- lapply(imgs, image_read)
img_joined <- image_join(img_list)
img_animated <- image_animate(img_joined, fps = 20)
#export en gif
image_write(image = img_animated, path = "Smoothed_binned_colored_animated_cartogram_Export/animation.gif")
include_graphics(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram_Export/animation.gif"))
#export en mp4
rootwd <- getwd()
setwd(paste0(getwd(),"/Smoothed_binned_colored_animated_cartogram"))
av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")
## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Smoothed_binned_colored_animated_cartogram\\animation.mp4"
file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Smoothed_binned_colored_animated_cartogram_Export"))
dir.create("Smoothed_continuous_colored_animated_cartogram_Export")
#chargement des informations des fichiers à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"), full.names = TRUE))
#tri par date de création
details = details[with(details, order(as.POSIXct(ctime))), ]
list <- rep("0",110)
for (i in 1:10){
list[i] <- paste("00",i, sep = "")
}
for (i in 10:110){
if (i>=10 & i<100){
list[i] <- paste("0",i, sep = "")
} else {
list[i] <- paste(i)
}
}
file.rename(rownames(details), paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram/anim", list[1:110],".png"))
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE
# rechargement des informations des fichiers renommés à intégrer à l'animation
details = file.info(list.files(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"), full.names = TRUE))
details <- details[15:110,]
imgs = rownames(details)
img_list <- lapply(imgs, image_read)
img_joined <- image_join(img_list)
img_animated <- image_animate(img_joined, fps = 20)
#export en gif
image_write(image = img_animated, path = "Smoothed_continuous_colored_animated_cartogram_Export/animation.gif")
include_graphics(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram_Export/animation.gif"))
#export en mp4
rootwd <- getwd()
setwd(paste0(getwd(),"/Smoothed_continuous_colored_animated_cartogram"))
av_encode_video(input = list.files(getwd()),framerate = 15, "animation.mp4")
## [1] "F:\\Rennes\\Recherche\\ANR MoDurAL\\Valorizacion\\IJC\\version anglaise\\Derniers_traitements\\Envoi Florent 12 01\\Smoothed_continuous_colored_animated_cartogram\\animation.mp4"
file.move(paste0(getwd(),"/animation.mp4"),paste0(rootwd,"/Smoothed_continuous_colored_animated_cartogram_Export"))